clc
clear all

Inpath='/Volumes/DPR/Data/IMERG/ERE_huffman/monthly/uncal_correct_date/';
Indir=dir(fullfile(Inpath,'*.mat'));
Inlength=length(Indir);
Outpath='/Volumes/DPR/Data/IMERG/ERE_huffman/threshold/threshold_all_uncal.mat';

temp_sum = zeros(242,602);
temp_count = zeros(242,602);

%overall mean
for i=1:Inlength
    i
    filename = strcat(Inpath,Indir(i).name);
    load(filename)
    
    meanRainfall(isnan(meanRainfall))=0;
    
    temp_count = temp_count + countRainfall;
    temp_sum = temp_sum + countRainfall.*meanRainfall;    
end

overall_mean = temp_sum./temp_count;

b = temp_count .* overall_mean .* overall_mean;

a = zeros(242,602);

for i=1:Inlength
    i
    filename = fullfile(Inpath,Indir(i).name);
    load(filename)
    a = a + sumOfSquareRainfall - 2 * overall_mean .* countRainfall .* meanRainfall;
end

variance = (a + b) ./ (temp_count - 1);
overall_std = sqrt(variance);

ERE_threshold = overall_mean + 2.5 * overall_std;

ERE_threshold(ERE_threshold<5) = 5;

ERE_threshold(isnan(ERE_threshold))=5;

save(Outpath,'ERE_threshold');

ERE_threshold = flipud(ERE_threshold);

percentage_raintime = (temp_count/70128);

figure
imagesc(ERE_threshold);
caxis([5,12]);
title('ERE Precip Threshold Uncalibrated IMERG 2014.4 to 2018.3','fontsize',12,'fontWeight','normal');
set(gca,'yticklabel',{''},'fontsize',8);
set(gca,'xticklabel',{''},'fontsize',15);
colorbar;

figure
imagesc(flipud(overall_mean));
title('ERE Precip overall mean IMERG 2018','fontsize',15,'fontWeight','normal');
set(gca,'yticklabel',{''},'fontsize',8);
set(gca,'xticklabel',{''},'fontsize',15);
%caxis([5,11]);
colorbar;

figure
imagesc(flipud(percentage_raintime));
title('ERE Percentage of Rain Time','fontsize',15,'fontWeight','normal');
set(gca,'yticklabel',{''},'fontsize',8);
set(gca,'xticklabel',{''},'fontsize',15);
caxis([0,0.25])
colorbar